Teoria Macroeconômica II - Artigo parte I

Dúvidas

1. Tem como fazer alguma conversão dos valores pré-plano real para tributo do trabalho e capital para não perder as 30 primeiras observações?

Aisha: Não tem o que fazer.

2. Tem que fazer alguma coisa com a sazonalidade e tendência?

Aisha: É para usar o filtro automático mas explicar bem o que ele faz.

3. No artigo de bayesiana eu havia utilizado a SELIC overnight. Aí fui ver no site do Bacen e só tem a série 4390- Taxa de juros (Selic acumulada no mês). Qual a diferença entre essas duas?

Aisha: Usar a SELIC do Bacen é melhor (depois tem que voltar aqui e escrever quais as diferenças entre essas séries e porque a do Bacen é melhor), porém tem que lembrar de passar para valores anuais (por enquanto estava mensal). Update: Fui procurar e entendi o que é a Selic Overnight (que tem no IPEA), mas ainda não entendi o que é a Selic acumulada no mês do Banco Central. Update: A Selic Acumulada no site do Bacen e a Selic Overnight do site do IPEA tem os mesmos valores. E a Selic Acumulada é feita usando a Selic diária (série 11) com a fórmula de juros compostos. Agora resta saber por que a série 11 tem formato tão de escadinha.

4. Quais outras variáveis adicionar no modelo?

Aisha: Olhar o que o Mumtaz usou no VAR dele. Basicamente procurar por produto, inflação, câmbio (talvez)…

5. Qual a temporalidade dos dados?

Aisha: De acordo com o Guilherme, se usar câmbio o melhor é fazer a partir do fim da âncora cambial. Fui olhar e isso implica começar a série em Janeiro/99, quando foi adotado o regime de metas de inflação (http://www.rep.org.br/pdf/87-1.pdf). Por enquanto eu vou rodar desde 94 e depois eu adapto, se for o caso. Fui ver e a série do IBC-Br começa em 2003, daí fode tudo.

6. O que é o equivalente tupiniquim de three month treasury bill rate?

Aisha: Dá para usar isso? https://cran.r-project.org/web/packages/GetTDData/vignettes/gtdd-vignette_GetTDData.html . Update: O Portela disse o seguinte " Oi Aisha, Você poderia usar a taxa do swap DI-Pré de 3 meses. Dá pra baixar do banco central ou então pelo BETS. Um abraço“. No banco central, a série é 7818- Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período). Mas eu não entendi que trem é esse… como que isso se relaciona com política monetária? E qual seria a justificativa para usar isso e não a taxa de juros da economia para colocar no VAR do Mumtaz?

7. Log differences é a diferença dos logs ou o log da diferença?

Aisha: Dar uma olhada no material de macro, uma delas era taxa de crescimento (diferença dos logs).

8. Como calcular a taxa de crescimento da taxa de câmbio efetiva nominal?

Aisha: Sei lá. O Mumtaz usa growth of the nominal effective exchange rate, mas aí não sei se ele faz a diferença dos logs e por isso chama assim ou se ele calcula a taxa de crescimento e depois faz a diferença dos logs. Eu posso dar uma olhada no código deles depois.

Baixando os dados

Dependências dos pacotes

library(ggplot2, quietly = TRUE)
library(forecast, quietly = TRUE)
library(BETS, quietly = TRUE)
library(seasonal, quietly = TRUE)
library(seasonalview, quietly = TRUE)
library(lubridate, quietly = TRUE)
library(zoo, quietly = TRUE)
library(stargazer, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(reshape2, quietly = TRUE)
library(ggfortify, quietly = TRUE)
library(RColorBrewer, quietly = TRUE)
library(scales, quietly = TRUE)
library(quantmod, quietly = TRUE)
library(PerformanceAnalytics, quietly = TRUE)
library(strucchange, quietly = TRUE)

Aisha: Caso o BETS não esteja instalado, retirar o eval=FALSE do chunk abaixo e rodar ele.

library(devtools)
devtools::install_github("pedrocostaferreira/BETS", force= TRUE)

Download dos dados

Os dados podem ser pegos de https://www3.bcb.gov.br/sgspub/localizarseries/localizarSeries.do?method=prepararTelaLocalizarSeries:

Renda do trabalho: Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (7620)

Renda do capital: Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (7621)

Ambas séries são mensais e começam em 31/01/1992. O último dado disponível em 05/01/2018 era de novembro de 2017, totalizando 311 meses. O pacote BETS consegue fazer download das séries apenas usando o código delas. Cortando em 1994 (por causa do plano real) ficam 281 observações (vou deixar para ver o trem da âncora cambial depois).

Update: Por enquanto estou incluindo IBC-Br o que me deixa com 177 observações apenas…

trabalho <- BETS.get("7620")
capital <- BETS.get("7620")

autoplot(trabalho)

autoplot(capital)

trabalho <- BETS.get("7620", from = "1994-07-01", to = "2017-11-01")
autoplot(trabalho)

capital <- BETS.get("7621", from = "1994-07-01", to = "2017-11-01")
autoplot(capital)

capital_trabalho <- capital/trabalho
autoplot(capital/trabalho)

length(capital_trabalho)

## [1] 281

Inicialmente eu fiz um primeiro VAR apenas com Selic e razão capital/trabalho. Havia a dúvida se era melhor usar a Selic do site do BACEN ou a Selic Overnight do site do IPEA. De acordo com o Guilherme, o melhor é o do Bacen pois é a efetiva (não entendi bem essa história…), mas lá em Bayesiana ele não me disse isso “porque estava ocupado me ensinando o que era Selic”.

Uma digressão sobre as taxas de juros SELIC

Do site do Banco Central, tem essa definição: * “Define-se Taxa Selic como a taxa média ajustada dos financiamentos diários apurados no Sistema Especial de Liquidação e de Custódia (Selic) para títulos federais. Para fins de cálculo da taxa, são considerados os financiamentos diários relativos às operações registradas e liquidadas no próprio Selic e em sistemas operados por câmaras ou prestadores de serviços de compensação e de liquidação” (Link).

Aqui tem uma explicação do que é a Selic Overnight:

• “(…) A Selic Over é obtida a partir do financiamento no mercado interbancário lastreado em títulos públicos, sendo calculada diariamente. (…) Essa taxa surge da necessidade de financiamento de quem está comprando títulos públicos.” No site ainda tem um exemplo que ajuda a entender melhor.

Mas isso ainda não explica o que é a série 4390 - Taxa de juros - Selic acumulada no mês que tem no sistema do Banco Central. Olhando na caixinha de explicações metodológicas, diz que ela é derivada da série 11 - Taxa de Juros - Selic. Essa taxa 11 está em %a.d.. Então vou pegar uns valores para ver o que acontece:

selic_11   <- BETS.get("11",   from = "2000-01-01", to = "2017-12-31")
#head(selic_11)
selic_4390 <- BETS.get("4390", from = "2000-01-01", to = "2017-12-31")
#head(selic_4390)

df1 <- as.data.frame(selic_11)
#head(df1)

p1 <- ggplot(df1, aes(date, value)) +
        geom_line() +
        scale_x_date() +
        xlab("Data") +
        ylab("Selic (%a.d.)") +
        ggtitle("11 - Taxa de Juros - Selic")


df2 <- data.frame(time(selic_4390), as.numeric(selic_4390))
#head(df2)
names(df2) <- c("date", "value")

p2 <- ggplot(df2, aes(date, value)) +
        geom_line() +
        xlab("Data") +
        ylab("Selic (%a.m.)") +
        ggtitle("4390 - Taxa de Juros - Selic acumulada no mês")

grid.arrange(p1, p2, ncol=1, nrow = 2)

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

Depois disso, eu peguei ambas séries de janeiro de 2017 a dezembro de 2017 e na série diária eu somei os valores de cada mês, para obter 12 valores e comparei com os 12 valores da série 4390. Os valores são bem similares. Depois eu tentei fazer um puxadinho do que é a fórmula nos dados básicos da série 4390, onde diz CONVPER((ACMVALORES(((SERIE(11)/100)+1),"mensal","multiplicacao")-1)*100,"mensal","ultimovalor"). O que eu fiz no Excel foi pegar cada valor, dividir por \(100\) e somar \(1\). Depois disso, multipliquei os valores dentro de um mesmo mês, subtraí \(1\) e multipliquei por \(100\).

Comparação da série 11 com a série 4390

Resumindo, o que tem na série 4390 parece ser a fórmula de juros compostos usando os valores da série 11, que por sua vez tem valores do tipo “escadinha”.

Agora, vou olhar para a Selic Overnight, do site do IPEA.

selic_overnight <- read.csv2("C:\\Users\\Aishameriane\\OneDrive\\Documentos\\Mestrado Economia\\Teoria Macroeconômica II - 2017-2\\Artigo\\Aplicação\\selic_overnight.csv", sep = ";", dec = ",", header = TRUE)
selic_overnight <- selic_overnight[,-3]
colnames(selic_overnight) <- c("Data", "Selic")

# Arrumando as datas
selic_overnight[,1]<-paste(selic_overnight[,1], ".01", sep="")
selic_overnight[,1]<-ymd(selic_overnight[,1])

inicio <- which(selic_overnight[,1] == "1994-07-01")
fim <- which(selic_overnight[,1] == "2017-11-01")
selic_overnight<-selic_overnight[inicio:fim,]

p3 <- ggplot(selic_overnight, aes(Data, Selic)) +
        geom_line() +
        xlab("Data") +
        ylab("Selic (%a.m.)") +
        ggtitle("Taxa de Juros - Overnight Selic")

grid.arrange(p1, p2, p3, ncol=1, nrow = 3)

Eu estou desconfiada que a Overnight e a série 4390 são a mesma coisa, então vou dropar uns valores da overnight e ver o que acontece.

inicio <- which(selic_overnight[,1] == "2000-01-01")
fim <- which(selic_overnight[,1] == "2017-11-01")
selic_overnight_menor<-selic_overnight[inicio:fim,]

p4 <- ggplot(selic_overnight_menor, aes(Data, Selic)) +
        geom_line() +
        xlab("Data") +
        ylab("Selic (%a.m.)") +
        ggtitle("Taxa de Juros - Overnight Selic")

grid.arrange(p1, p2, p4, ncol=1, nrow = 3)

Os gráficos são bem similares… Mas e os valores? (obs: eu olhei para os valores da cauda também e usei janelas maiores de dados… só os dois últimos não são iguais na segunda casa decimal.)

cbind(head(round(selic_overnight_menor[,2],2), 10), head(selic_4390[1:(length(selic_4390)-1)], 10))

##       [,1] [,2]
##  [1,] 1.46 1.46
##  [2,] 1.45 1.45
##  [3,] 1.45 1.45
##  [4,] 1.30 1.30
##  [5,] 1.49 1.49
##  [6,] 1.39 1.39
##  [7,] 1.31 1.31
##  [8,] 1.41 1.41
##  [9,] 1.22 1.22
## [10,] 1.29 1.29

E acho que isso resolve uma parte do mistério. A série 4390 e a série Selic Overnight do site do IPEA de fato parecem ser a mesma série. E elas são o resultado da fórmula de juros compostos na série 11. Aisha: Agora resta saber por que a série 11 tem formato tão de escadinha.

Como é mais fácil de pegar os dados usando o pacote BETS, eu vou usar a série selic_4390 como a série “oficial”. Agora resta baixar os dados da série inteira (jul-94 até dez-17) e anualizar.

selic_4390 <- BETS.get("4390", from = "1994-07-01", to = "2017-12-31") 

## Transformando em série anual
selic <- ((1+selic_4390/100)^(12)-1)*100
head(selic)

##            Jul       Aug       Sep       Oct       Nov       Dec
## 1994 121.95744  63.27210  56.99082  53.22268  61.40116  56.44736

tail(selic)

##            Jul       Aug       Sep       Oct       Nov       Dec
## 2017 10.033869 10.033869  7.956187  7.956187  7.058561  6.675963

length(selic)

## [1] 282

df2 <- data.frame(time(selic), as.numeric(selic))
names(df2) <- c("Data", "Selic")

p5 <- ggplot(df2, aes(Data, Selic)) +
        geom_line() +
        xlab("Data") +
        ylab("Selic (%a.a.)") +
        ggtitle("Taxa de Juros - Selic")

p5

Os valores do início da série estão meio ruins, mas depois eu me preocupo com isso se for o caso. Tem outra coisa me incomodando que os valores estão meio diferentes do que tem no site do Bacen: https://www.bcb.gov.br/Pec/Copom/Port/taxaSelic.asp. Aí não sei se eu não fiz merda na hora de converter a taxa.

Fim da digressão sobre a SELIC

Pegando os outros dados

As variáveis do modelo do Mumtaz são:

• PIB per capita (real GDP per capita) - peguei IBC-Br 24364 que já tem ajuste sazonal (problema é que a série começa em janeiro de 2003 e acaba em outubro de 2017). Depois acabei pegando o PIB Mensal do Ipea Data deflacionado usando o IPCA de novembro de 2017. Essa deflação é feita com alguma bruxaria que não consegui identificar e por isso estou carregando o arquivo csv - não deu para baixar as séries do Bacen e deflacionar eu mesma e achar a mesma coisa.
• Inflação (CPI inflation) - peguei IPCA 433 que já está em variação percentual mensal;
• Three month treasury bill rate - falei com o Portela e ele sugeriu usar a série 7818 Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) (ela tem valores de 30/09/1999 a desembro de 2017). Depois comparei com a selic e como são praticamente iguais, fiquei com a Selic, mas fiz um teste com ela.
• Taxa de crescimento da taxa de câmbio efetiva nominal (growth of the nominal effective exchange rate) - mínima ideia de como calcular isso também. Peguei a série que o Portela usava em séries temporais: 3696 - taxa de câmbio livre - dólar americano (venda) fim de período mensal.

Todas as variáveis exceto as de desigualdade, inflação e a taxa de juros estão em log differences (diferença dos logaritmos).

PIB Per capita

Bom, não tem PIB per capita, então tem que ser IBC-Br.

ibcbr <- BETS.get("24364", from = "2003-09-01", to = "2017-12-31") 

df3 <- data.frame(time(ibcbr), as.numeric(ibcbr))
names(df3) <- c("Data", "IBCBr")

p6 <- ggplot(df3, aes(Data, IBCBr)) +
        geom_line() +
        xlab("Data") +
        ylab("IBC-Br") +
        ggtitle("Índice de Atividade Econômica do Bacen - com ajuste sazonal")

p6

LEMBRETE: Esse trem tem uma tendência monstra, depois tem que ver isso. Tirei a diferença da série em log.

LEMBRETE: A série começa em janeiro de 2003… Daí fica com 178 dados… tem que ver com o Guilherme depois.

PIB Mensal

Fuçando no site do Banco Central, eu achei uma série chamada PIB mensal - Valores correntes (R$ milhões) (4380).

Vou dar uma olhada nela:

PIBmensal <- BETS.get("4380", from = "1994-07-01", to = "2017-12-31") 
data1 <- seq(as.Date("1994-07-01"), length = length(PIBmensal), by = "1 month")

df3 <- data.frame(data1, as.numeric(PIBmensal))

p6 <- ggplot(df3, aes(data1, PIBmensal)) +
        geom_line() +
        xlab("Data") +
        ylab("PIB Mensal (R$ milhões)") +
        ggtitle("PIB mensal - Valores correntes")

p6

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

# Como eu estou com preguiça de calcular o PIB deflacionado, peguei a série do site do IPEA pronta deflacionada pelo IPCA bonitinha

PIBmensal2 <- read.csv("C:\\Users\\Aishameriane\\OneDrive\\Documentos\\Mestrado Economia\\Teoria Macroeconômica II - 2017-2\\Artigo\\Aplicação\\PIBmensal.csv", header = T, sep = ";", dec = ",")
data1 <- seq(as.Date("1994-07-01"), length = length(PIBmensal2), by = "1 month")

df3 <- data.frame(data1, as.numeric(PIBmensal))
names(df3) <- c("Data", "PIBmensal")

p6 <- ggplot(df3, aes(Data, PIBmensal)) +
        geom_line() +
        xlab("Data") +
        ylab("PIB Mensal (R$ milhões)") +
        ggtitle("PIB mensal - Valores reais (base IPCA 11/2017)")

p6 + scale_y_continuous(labels = comma)

# Agora preciso transformar em índice
PIBmensal_in <- (PIBmensal/PIBmensal[1])*100
df3 <- data.frame(data1, as.numeric(PIBmensal_in))
names(df3) <- c("Data", "PIBmensal")

p6 <- ggplot(df3, aes(Data, PIBmensal)) +
        geom_line() +
        xlab("Data") +
        ylab("PIB Mensal") +
        ggtitle("PIB mensal índice - base = 2004/07")

p6

ibcbr <- BETS.get("24364", from = "2003-01-01", to = "2017-10-31")
data2 <- seq(as.Date("2003-01-01"), length = length(ibcbr), by = "1 month")

df4 <- data.frame(data2, as.numeric(ibcbr))
names(df4) <- c("Data", "IBCBr")

# Tentando juntar os pontos
test1 <- melt(data = df3, id.vars = "Data")
test2 <- melt(data = df4, id.vars = "Data")

test <- rbind(test1, test2)

ggplot(test, aes(Data, value, colour = variable)) +
  geom_line(alpha = 1)+
  labs(title="IBC-Br x PIB mensal", y = "", x= "Tempo", color = "Variável") +
  scale_colour_brewer(palette = "Dark2") +
  theme_bw()

PIBmensal_in <- ts(PIBmensal_in,  start = c(1994, 07, 01), frequency = 12) 

Update: Eu fiz alguns testes e o VAR com o PIB mensal “deu bom”. Como é importante saber fazer o que um aluno graduado em economia faz (MOURA, 2018), vou refazer as coisas em cima deflacionando eu mesma a série. Peguei um tutorial aqui. Update: Essa merda não funciona, vou continuar na ignorância mesmo, aff. Update: Havia um erro na merda da série do IPCA do IPEA data e aí os resultados entre o que eu calculava no Excel e o que tinha no Bacen davam diferentes. Eu e o Cássio mandamos um e-mail pro IPEA avisando. E agora eu consigo fazer tudo sem usar um csv. :D

Pra Aisha do futuro, eu tentei usar a fórmula e os valores deram diferentes do site do IPEA. Aí lendo as notas metodológicas, fala isso aqui: " Os índices IPCA, INPC e IGP-DI refletem o nível de preços médios de cada período, sendo inadequados para deflacionar séries referentes ao final do período. Visando minimizar tal problema, o sistema calcula esses índices de preço no final de cada período através do cálculo da média geométrica entre o índice do período e o índice do período seguinte." Calcular a média geométrica é tirar a raiz quadrada do produto, mas quando o índice dá negativo, dá ruim. Aí tentei fazer a média geométrica do IPCA acumulado, também não deu e eu desisti porque não tá sobrando tanto tempo assim. Agora deu.

PIBmensal <- BETS.get("4380", from = "1999-06-01", to = "2017-11-30") 
data1 <- seq(as.Date("1994-07-01"), length = length(PIBmensal), by = "1 month")
df3 <- data.frame(data1, as.numeric(PIBmensal))
names(df3) <- c("Data", "PIBmensal")

p6 <- ggplot(df3, aes(Data, PIBmensal)) +
        geom_line() +
        xlab("Data") +
        ylab("PIB Mensal (R$ milhões)") +
        ggtitle("PIB mensal - Valores correntes")
p6

# O IPCA tem que ser pego desde 1993 caso queira confrontar com os valores do IPEA Data
ipca <- BETS.get("433", from = "1993-12-31", to = "2017-12-31") 
ipca2 <- ipca
ipca2[1] <- 100

# Transformar em índice
for (i in 2:length(ipca2)){
  ipca2[i] <- round(ipca2[(i-1)]*(1+ipca2[(i)]/100),2)
}

# Média geométrica
for (i in 1:(length(ipca2)-1)){
  ipca2[i] <- sqrt(ipca2[i]*ipca2[(i+1)])
}

ipca2 <- ipca2[1:(length(ipca2)-1)]

# Mudança para base de nov/2017
for (i in 1:length(ipca2)) {
  ipca2[i] <- tail(ipca2, n=1)/ipca2[i]
}

# Tirando as observações que não quero
ipca2 <- ts(ipca2,  start = c(1993, 12, 31), frequency = 12)
ipca2 <- window(ipca2, c(1996, 6))

PIBmensal <- PIBmensal*ipca2

Inflação

Para inflação usei a série 433, do IPCA. Tem duas coisas que podem ser feitas, uma é calcular o índice acumulado no ano (o que não parece interessante porque implica que todo início de ano teria uma queda) e a segunda é calcular o índice acumulado nos últimos 12 meses. Para fazer a segunda, tem que primeiro pegar todos os valores da série e dividir por 100 e depois acrescentar 1. Dessa nova série, são multiplicados os 12 valores anteriores (incluindo o atual), desconta um e multiplica por 100.

# Primeiro um gráfico simples
ipca <- BETS.get("433", from = "1994-07-01", to = "2017-12-31") 

df4 <- data.frame(time(ipca), as.numeric(ipca))
names(df4) <- c("Data", "IPCA")

p7 <- ggplot(df4, aes(Data, IPCA)) +
        geom_line() +
        xlab("Data") +
        ylab("IPCA") +
        ggtitle("Inflação (IPCA) - mensal")

p7

# Agora a inflação acumulada em 12 meses
ipca_acum <- ipca/100 + 1
ipca_acum2 <- vector()

for (i in 12:282){
  ipca_acum2[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}
ipca_acum2 <- ts(ipca_acum2,  start = c(1995, 6, 1), frequency = 12)

df4 <- data.frame(time(ipca_acum2), as.numeric(ipca_acum2))
names(df4) <- c("Data", "IPCA")

p8 <- ggplot(df4, aes(Data, IPCA)) +
        geom_line() +
        xlab("Data") +
        ylab("IPCA") +
        ggtitle("Inflação (IPCA) - acumulada 12 meses")

p8

Eu conferi os valores com isso aqui e bateu. Já posso pedir a carteirinha de economista, yay! Posso nada, não consegui deflacionar o PIB. CONSEGUI, mwahahah.

Aisha: De novo tem aquelas valores ruins no início da série, acho que vou ter mesmo que eliminar uns valores.

Three month treasury bill rate

Depois de falar com o Portela, cheguei nessa série aqui: 7818 Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) . Eu ainda não sei como que isso se relaciona a uma taxa de juros (no sentido de instrumento de política monetária), mas vou deixar os dados aqui guardados por enquanto.

swap90 <- BETS.get("7818", from = "1999-09-01", to = "2017-12-31")

df7 <- data.frame(time(swap90), as.numeric(swap90))
names(df7) <- c("Data", "DI90")

p10 <- ggplot(df7, aes(Data, DI90)) +
        geom_line() +
        xlab("Data") +
        ylab("DI90 (%a.a.)") +
        ggtitle("Taxa referencial de swaps DI pré 3 meses - fim de período")

p10

Vou tentar comparar com a SELIC. Será que essas duas taxas são comparáveis? Acho que sim porque eu passei SELIC para %a.a. também….

selic_4390 <- BETS.get("4390", from = "1999-09-01", to = "2017-12-31") 

## Transformando em série anual
selic <- ((1+selic_4390/100)^(12)-1)*100
df10 <- data.frame(time(selic), as.numeric(selic))
names(df10) <- c("Data","Selic")
p11 <- ggplot(df10, aes(Data, Selic)) +
        geom_line() +
        xlab("Data") +
        ylab("Selic (%a.a.)") +
        ggtitle("Selic")

grid.arrange(p10, p11, ncol=1, nrow = 2)

df8 <- data.frame(df7, selic)
names(df8) <- c("Data", "DI90", "Selic")

df9 <- melt(data = df8, id.vars = "Data")

ggplot(df9, aes(Data, value, colour = variable)) +
  geom_line(alpha = 1)+
  labs(title="Taxas de juros", y = "", x= "Time", color = "Taxa") +
  scale_colour_brewer(palette = "Dark2") +
  theme_bw()

Elas são praticamente iguais, então vou ficar com a Selic. Update (18/01/2018) Guilherme mandou ficar com o swap DI. Update (22/01/2018) Agora que vi que essa coisa começa só em 99.

Câmbio

Do câmbio, tem essas séries no Bacen:

• Taxa de câmbio - Livre - Dólar americano (compra) - Média de período - mensal: número da série 3697, peridiocidade mensal, unidade u.m.c./US$.
• Taxa de câmbio - Livre - Dólar americano (venda) - Fim de período - mensal: número da série 3696, peridiocidade mensal, unidade u.m.c./US$.

Como o Portela usa a segunda, vou usar a segunda (argumento de autoridade).

cambio <- BETS.get("3696", from = "1994-07-01", to = "2017-12-31")

df5 <- data.frame(time(cambio), as.numeric(cambio))
names(df5) <- c("Data", "Cambio")

p8 <- ggplot(df5, aes(Data, Cambio)) +
        geom_line() +
        xlab("Data") +
        ylab("Tx. de Câmbio (u.m.c./US$)") +
        ggtitle("Taxa de câmbio - Dólar americano (venda) - Fim de período - mensal")

p8

Realmente a partir de 99 que as coisas começaram a ficar diferentes.

Vou tentar calcular a taxa de crescimento usando a diferença dos logaritmos.

cambio_cres <- diff(log(cambio), 1)

df6 <- data.frame(time(cambio_cres), as.numeric(cambio_cres))
names(df6) <- c("Data", "Cambio")

p9 <- ggplot(df6, aes(Data, Cambio)) +
        geom_line() +
        xlab("Data") +
        ylab("Variação") +
        ggtitle("Taxa de variação da Taxa de câmbio - Dólar americano (venda) - Fim de período - mensal")

p9

Ibovespa

Aí como estava com tempo sobrando, pensei em pegar o Ibovespa. Teoricamente teria alguma relação entre preço dos ativos e taxa de juros (esse é um dos canais de desigualdade).

symbols <- c("^BVSP")
getSymbols(symbols, src = 'yahoo', from = "1996-01-01", 
             auto.assign = TRUE, warnings = FALSE) 

## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.

## 
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
## 
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).

## Warning: ^BVSP contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.

## [1] "BVSP"

prices <- merge(BVSP)
prices <- prices[,grepl( "Close" , names(prices) )]
colnames(prices) <- symbols
returns <- na.omit(monthlyReturn(prices, type = "log"))

## Warning in to_period(xx, period = on.opts[[period]], ...): missing values
## removed from data

plot.xts(returns)

ibovespa <- ts(returns, start = c(1996,1,1), frequency = 12)

df1 <- data.frame(seq(as.Date("1996-01-01"), length = length(ibovespa), by = "1 month"), ibovespa)
names(df1) <- c("Data", "Ibovespa")

p9 <- ggplot(df1, aes(Data, Ibovespa)) +
        geom_line() +
        xlab("Data") +
        ylab("Ibovespa") +
        scale_x_date(date_breaks = "12 months")+
        ggtitle("Log retornos mensais do Ibovespa")

p9 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

Todos os dados juntos

A lista de dados que tenho até agora é:

1. Razão capital-trabalho (para desigualdade), calculada a partir das seguintes séries:
   • Renda do trabalho: Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (7620)
   • Renda do capital: Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (7621) Não foi feita nenhuma transformação nos dados.
2. Taxa Selic (em substituição à three month treasury bill rate) calculada a partir da série:

   • Taxa de juros - Selic acumulada no mês (4390) Eu comparei com a taxa referencial de swaps DI pré-fixada e ambas são praticamente a mesma coisa, então fiquei com a Selic mesmo. A única transformação feita (até agora) foi passar de taxa mensal para anual usando a fórmula: \(\left((1+tx/100)^12 -1\right)*100\).

   • Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) (7818) Sem transformação (a taxa já é ao ano). O problema é que essa série começa apenas em 1999.

3. IPCA (para inflação), utilizando a seguinte série:
   • Índice nacional de preços ao consumidor-amplo (IPCA) (433) Como ele está em variação percentual mensal, foi necessário transformar para acumulado dos últimos 12 meses utilizando a fórmula: \[IPCA_i = \left[\left(\prod\limits_{j=i-11}^i \left(\frac{IPCA_j}{100}+1\right) \right) -1\right]*100\]
4. Para o PIB per capita
   • IBC-Br (para o PIB per capita), utilizando a seguinte série:
     • Índice de Atividade Econômica do Banco Central (IBC-Br) - com ajuste sazonal (24364) Não fiz nenhuma transformação por enquanto, mas o bicho tem uma baita tendência. Fiz a diferença da série em log, igual o câmbio.
   • PIB mensal (retirado do site do IPEA) deflacionado usando o IPCA de novembro de 2017. Fiz a diferença da série em log, igual o câmbio.
5. Crescimento da taxa de câmbio (para growth of the nominal effective exchange rate). Usei a série:
   • Taxa de câmbio - Livre - Dólar americano (venda) - Fim de período - mensal (3696) Como o Mumtaz disse que usou log diferença, eu calculei o log e tirei primeira diferença nessa série. As outras eu ainda não fiz, porque o IBC-Br é índice e a inflação também…

rm(list = ls())

# Auxiliary variables, so I don't need to bother when something changes
inicio <- "2003-02-01"
fim    <- "2017-10-31"

# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca   <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries

trabalho <- BETS.get("7620", from = inicio, to = fim)
capital  <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho

# Usando ndiffs(ibcbr,test="adf",alpha = 0.1) se encontra que o ibcbr é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
ibcbr_raw    <- BETS.get("24364", from = inicio_cambio, to = fim)
ibcbr        <- diff(log(ibcbr_raw), 1)

cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio     <- diff(log(cambio_raw), 1)

selic_4390 <- BETS.get("4390", from = inicio, to = fim) 
selic <- ((1+selic_4390/100)^(12)-1)*100

ipca_raw <- BETS.get("433", from = inicio_ipca, to = fim) 
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()

final <- length(ipca_acum)
for (i in 12:final){
  ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}

ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))

ipca <- ts(ipca,  start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD

# Plotting some nice graphs

df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho, selic, ibcbr, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "IBCBr", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")


ggplot(df2, aes(Data, value, colour = variable)) +
  geom_line(alpha = 1)+
  labs(title="Variáveis do modelo", y = "", x= "Time", color = "Variável") +
  scale_colour_brewer(palette = "Dark2") +
  theme_bw()

cores <- brewer.pal(5, "Dark2")

# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Capital trabalho") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p2 <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[2])+
        scale_y_continuous(name="Selic") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p3 <- ggplot(df2[which(df2$variable == "IBCBr"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[3])+
        scale_y_continuous(name="IBC-Br") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[4])+
        scale_y_continuous(name="Tx. Câmbio (var)") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p5 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[5])+
        scale_y_continuous(name="IPCA") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
        theme_bw()
p5 <- p5 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5)

Acrescentando o PIB e o Swap DI 3 meses:

rm(list = ls())

# Auxiliary variables, so I don't need to bother when something changes
inicio <- "1996-01-01"
fim    <- "2017-11-30"
inicio_deflator <- "1993-12-31"

# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca2   <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries
fim_ipca      <- paste(seq(as.Date(fim), length = 2, by = "1 month")[2]) # 1 mês depois do fim das outras séries

trabalho <- BETS.get("7620", from = inicio, to = fim)
capital  <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho

ipca_raw <- BETS.get("433", from = inicio_deflator, to = fim_ipca) 
ipca2 <- ipca_raw
ipca2[1] <- 100

# Transformar em índice
for (i in 2:length(ipca2)){
  ipca2[i] <- round(ipca2[(i-1)]*(1+ipca2[(i)]/100),2)
}

# Média geométrica
for (i in 1:(length(ipca2)-1)){
  ipca2[i] <- sqrt(ipca2[i]*ipca2[(i+1)])
}

# Mudança para base de nov/2017
for (i in 1:length(ipca2)) {
  ipca2[i] <- tail(ipca2, n=2)[1]/ipca2[i]
}

# Tirando as observações que não quero
ipca2 <- ts(ipca2,  start = c(1993, 12, 31), frequency = 12)
ano_ipca <- as.numeric(substr(inicio_cambio, start = 1, stop = 4))
mes_ipca <- as.numeric(substr(inicio_cambio, start = 6, stop = 7))
ipca2 <- window(ipca2, c(ano_ipca, mes_ipca))

ano_ipca <- as.numeric(substr(inicio_ipca2, start = 1, stop = 4))
mes_ipca <- as.numeric(substr(inicio_ipca2, start = 6, stop = 7))
ipca_raw <- window(ipca_raw, c(ano_ipca, mes_ipca))
# Tira dezembro 2017
ipca_raw <- ipca_raw[1:(length(ipca_raw)-1)]

ipca_acum <- ipca_raw/100 + 1
ipca <- vector()

final <- length(ipca_acum)
for (i in 12:final){
  ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}

ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))

ipca <- ts(ipca,  start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD

# Usando ndiffs(PIBmensal,test="adf",alpha = 0.05) se encontra que o PIBmensal é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
PIBmensal <- BETS.get("4380", from = inicio_cambio, to = fim)
PIBmensal <- PIBmensal * ipca2 
PIB <- diff(log(PIBmensal), 1)

cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio     <- diff(log(cambio_raw), 1)

swap90 <- BETS.get("7818", from = inicio, to = fim)
swap90 <- ts.union(swap90, ipca)
swap90 <- swap90[,1]

selic_4390 <- BETS.get("4390", from = inicio, to = fim) 
selic <- ((1+selic_4390/100)^(12)-1)*100

# Plotting some nice graphs

df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho, swap90, PIB, cambio, selic, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Swap90", "PIB", "Cambio", "selic", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")

ggplot(df2, aes(Data, value, colour = variable)) +
  geom_line(alpha = 1)+
  labs(title="Variáveis do modelo", y = "", x= "Time", color = "Variável") +
  scale_colour_brewer(palette = "Dark2") +
  theme_bw()

cores <- brewer.pal(6, "Dark2")

# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Capital trabalho") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p2 <- ggplot(df2[which(df2$variable == "Swap90"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[2])+
        scale_y_continuous(name="SwapDI90 (%a.a.)") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p3 <- ggplot(df2[which(df2$variable == "PIB"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[3])+
        scale_y_continuous(name="PIB mensal (tx.cresc)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
        theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p3 <- p3 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[4])+
        scale_y_continuous(name="Tx. Câmbio (var)") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p5 <- ggplot(df2[which(df2$variable == "selic"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[5])+
        scale_y_continuous(name="Selic (%a.a)") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p5 <- p5 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

p6 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[6])+
        scale_y_continuous(name="IPCA") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
        theme_bw()
p6 <- p6 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))

grid.arrange(p1, p2, p4, p5, p3, p6, ncol=2, nrow = 3)

grid.arrange(p1, p4, p3, ncol=1, nrow = 3)

grid.arrange(p2, p5, p6, ncol = 1, nrow = 3)

A sazonalidade

Agora começa outro problema para tratar os dados: para a razão capital trabalho original, tem 30 observações de antes do plano real que tem valores muito altos. Depois de eliminar elas, tem duas séries com tendência e sazonalidade. A razão entre elas não aparenta ter tendência (fiz um teste que indicou que a série precisaria de uma diferença para ficar estacionária), porém parece ter uma componente sazonal persistente.

Resolvi tentar duas coisas: a primeira foi fazer manualmente a remoção de tendência e sazonalidade e a segunda foi usar um filtro pronto. Meu filtro manual ficou uma merda, então vou ficar com o automático mesmo.

Procedimento manual

Vou deixar aqui e qualquer coisa apago no final. Os chunks não estão sendo rodados. Usei este tutorial - https://anomaly.io/seasonal-trend-decomposition-in-r/.

Removendo tendência

Não sei se precisa. Como os dados são mensais, vou usar uma janela móvel de tamanho 12.

# Checa número de diferenciações necessárias para tornar séries estacionárias
ndiffs(capital_trabalho,test="adf",alpha = 0.1)
# De acordo com o teste, seria necessário diferenciar uma vez a série para que ela fique estacionária

# Remove a tendência ajustando uma janela móvel
trend_capital_trabalho = ma(capital_trabalho, order = 12, centre = T)
plot(as.ts(capital_trabalho))
lines(trend_capital_trabalho)
plot(as.ts(trend_capital_trabalho))

df <- as.data.frame(cbind(time(trend_capital_trabalho),trend_capital_trabalho, capital_trabalho))
head(df)
names(df) <- c("Data", "Tendencia", "CapitalTrabalho")
ggplot(df, aes(Data)) +
  geom_line(aes(y = CapitalTrabalho, colour = "Capital/Trabalho")) +
  geom_line(aes(y = Tendencia, colour = "Tendência")) 

Removendo sazonalidade

m_capital_trabalho = t(matrix(data = capital_trabalho, nrow = 12))
seasonal_capital_trabalho = colMeans(m_capital_trabalho, na.rm = T)
plot(as.ts(rep(seasonal_capital_trabalho,12)))

Residual

Vamos calcular o “random noise” usando \(Resíduo = Série - Sazonalidade - Tendência\).

residuo = capital_trabalho - seasonal_capital_trabalho - trend_capital_trabalho
plot(as.ts(residuo))

Ainda assim está meio estranho… Vou tentar usar um filtro automático.

X-13ARIMA-SEATS

O X-13ARIMA-SEATS é um algoritmo de ajuste sazonal desenvolvido pelo United States Census Bureau. Sua implementação no R é feita pelo pacote seasonal. Um manual pode ser visto aqui: https://cran.r-project.org/web/packages/seasonal/vignettes/seas.pdf. Usei muito este documento aqui: file:///C:/Users/Aishameriane/Desktop/Ferreira_Mattos_2016.pdf e também o teste de quebra estrutural que tem aqui https://pythonandr.com/2016/11/08/endogenously-detecting-structural-breaks-in-a-time-series-implementation-in-r/ .

Update (18/01/2018) O Guilherme sugeriu passar um filtro na quebra da razão capital trabalho em 1999. Primeiro vou tirar a sazonalidade da série toda, ver como fica o gráfico, quebrar em dois períodos a série original, depois ajudar o filtro sazonal para cada uma das séries.

#######################
# Capital trabalho    #
#######################
cores <- brewer.pal(6, "Dark2")

# Fazendo um teste incluindo os dados desde 1996 (que depois quero fazer o VAR com isso)
trabalho <- BETS.get("7620", from = "1996-01-01", to = "2017-11-01")
capital <- BETS.get("7621", from = "1996-01-01", to = "2017-11-01")
capital_trabalho <- capital/trabalho

m <- seas(x = capital_trabalho, transform.function = "none", regression.aictest = NULL)
qs(m)

##                     qs   p-val
## qsori        285.61105 0.00000
## qsorievadj   378.83903 0.00000
## qsrsd          0.00000 1.00000
## qssadj        46.77278 0.00000
## qssadjevadj    0.00000 1.00000
## qsirr         48.12373 0.00000
## qsirrevadj     0.00000 1.00000
## qssori       127.09650 0.00000
## qssorievadj  146.78226 0.00000
## qssrsd         0.00000 1.00000
## qsssadj       11.73690 0.00283
## qsssadjevadj   0.00000 1.00000
## qssirr         9.09815 0.01058
## qssirrevadj    0.00000 1.00000

plot(m)

summary(m)

## 
## Call:
## seas(x = capital_trabalho, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## Constant        0.023146   0.005041   4.592 4.40e-06 ***
## AO1996.Jan     -0.657991   0.085330  -7.711 1.25e-14 ***
## LS1996.Jul     -0.215914   0.035320  -6.113 9.77e-10 ***
## LS1998.Jan      0.851747   0.060992  13.965  < 2e-16 ***
## LS1998.Feb     -0.456248   0.075658  -6.030 1.64e-09 ***
## LS1998.Apr     -0.343722   0.057834  -5.943 2.79e-09 ***
## LS1998.Jul      0.339913   0.036980   9.192  < 2e-16 ***
## AO1998.Aug      0.373609   0.060735   6.151 7.68e-10 ***
## AO1999.Feb      0.712392   0.065041  10.953  < 2e-16 ***
## AO1999.Mar      0.518619   0.065041   7.974 1.54e-15 ***
## LS1999.Dec     -0.237452   0.024535  -9.678  < 2e-16 ***
## AO2001.Jun      0.297887   0.059118   5.039 4.68e-07 ***
## AO2001.Sep      0.267972   0.059117   4.533 5.82e-06 ***
## AO2001.Oct      0.301303   0.065105   4.628 3.69e-06 ***
## AO2002.Oct      0.580093   0.065104   8.910  < 2e-16 ***
## AO2003.Feb      0.330971   0.059131   5.597 2.18e-08 ***
## AO2003.Dec     -0.271940   0.060961  -4.461 8.16e-06 ***
## LS2004.Mar     -0.196492   0.025369  -7.745 9.53e-15 ***
## LS2004.Nov     -0.182872   0.026266  -6.962 3.35e-12 ***
## AO2005.Jun      1.001053   0.072560  13.796  < 2e-16 ***
## AO2005.Dec      0.394815   0.059238   6.665 2.65e-11 ***
## AO2006.Jan      0.320392   0.065114   4.920 8.63e-07 ***
## AO2006.Jun      0.996463   0.081624  12.208  < 2e-16 ***
## AO2007.Jan      0.265836   0.065064   4.086 4.39e-05 ***
## AO2007.Jun      0.691402   0.080924   8.544  < 2e-16 ***
## AO2008.Jun      0.244229   0.070170   3.481  0.00050 ***
## LS2009.Aug     -0.099893   0.023161  -4.313 1.61e-05 ***
## AO2011.Oct      0.261397   0.059084   4.424 9.68e-06 ***
## AO2011.Dec      0.241221   0.059084   4.083 4.45e-05 ***
## AO2013.Jun     -0.411605   0.065042  -6.328 2.48e-10 ***
## AO2014.Jun     -0.335137   0.065042  -5.153 2.57e-07 ***
## MA-Seasonal-12  0.163729   0.063369   2.584  0.00977 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 0 0)(0 1 1)  Obs.: 263  Transform: none
## AICc: -495.2, BIC: -389.2  QS (no seasonality in final):46.77 ***
## Box-Ljung (no autocorr.): 38.53 * Shapiro (normality): 0.9772 ***
## Messages generated by X-13:
## Notes:
## - Unable to test AO1998.Jan due to regression matrix
##   singularity.

sliding <- series(m, "sfs")

## specs have been added to the model: slidingspans

sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100

## [1] 0

ggmonthplot(capital_trabalho) + 
  theme_bw() +
  scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
  labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "Nenhuma transformação")

# Gráfico por ano e mês
ggseasonplot(capital_trabalho, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
  labs(title = "Média mensal da razão capital-trabalho", x = "Mês", subtitle = "Nenhuma transformação")

# Usa a interface gráfica
#view(m)

# Salva a série final em uma nova variável
capital_trabalho2 <- final(m)
#plot(capital_trabalho2)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho2), by = "1 month"), capital_trabalho, capital_trabalho2)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho, colour = "Original")) +
  geom_line(aes(y = capital_trabalho2, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
  scale_y_continuous(name = "Capital/Trabalho")+
  labs(color="Variável") +
  theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1

# Dividindo a série em pedaços distintos: 
# pega os pontos de quebras
bp_ts <- breakpoints(capital_trabalho ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = capital_trabalho ~ 1)
## 
## Breakpoints at observation number:
##                          
## m = 1      97            
## m = 2   39 97            
## m = 3   39 97         221
## m = 4   39 97 138     221
## m = 5   39 97 138 185 224
## 
## Corresponding to breakdates:
##                                                
## m = 1           2004(1)                        
## m = 2   1999(3) 2004(1)                        
## m = 3   1999(3) 2004(1)                 2014(5)
## m = 4   1999(3) 2004(1) 2007(6)         2014(5)
## m = 5   1999(3) 2004(1) 2007(6) 2011(5) 2014(8)
## 
## Fit:
##                                              
## m   0      1      2      3      4      5     
## RSS  22.23  20.07  19.79  19.63  19.42  19.40
## BIC 107.76  92.01  99.40 108.37 116.75 127.56

# Calcula o intervalo de confiança
ci_ts <- confint(bp_ts)
# Monta o dataframe
Data <- seq(as.Date("1996-01-01"), length = length(capital_trabalho), by = "1 month")
df1 <- melt(data.frame(Data,capital_trabalho), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "capital_trabalho"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Razão Capital/Trabalho") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(capital_trabalho), xend = Data[ci_ts$confint[3]], yend = min(capital_trabalho)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1

# Salvando duas séries para passar o seasonal
capital_trabalho_a <- capital_trabalho[1:bp_ts$breakpoints]
capital_trabalho_b <- capital_trabalho[(bp_ts$breakpoints+1):length(capital_trabalho)]
capital_trabalho_a <- ts(capital_trabalho_a,  start = c(1996, 01, 01), frequency = 12)
capital_trabalho_b <- ts(capital_trabalho_b,  start = c(2004, 02, 01), frequency = 12)

m <- seas(x = capital_trabalho_a, transform.function = "none", regression.aictest = NULL)
qs(m)

##                    qs   p-val
## qsori        34.12733 0.00000
## qsorievadj   91.53477 0.00000
## qsrsd         0.00000 1.00000
## qssadj        0.00000 1.00000
## qssadjevadj   0.00000 1.00000
## qsirr         0.00078 0.99961
## qsirrevadj    0.00000 1.00000
## qssori       33.13775 0.00000
## qssorievadj  85.86391 0.00000
## qssrsd        0.00000 1.00000
## qsssadj       0.00000 1.00000
## qsssadjevadj  0.00000 1.00000
## qssirr        0.05069 0.97497
## qssirrevadj   0.00000 1.00000

plot(m)

summary(m)

## 
## Call:
## seas(x = capital_trabalho_a, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AO1996.Jan        -0.35561    0.10377  -3.427 0.000611 ***
## AO1998.Jan         0.67046    0.08634   7.766 8.12e-15 ***
## LS1998.Jul         0.28775    0.05339   5.390 7.06e-08 ***
## AO1998.Aug         0.29547    0.08564   3.450 0.000561 ***
## AO1999.Feb         0.68086    0.09704   7.016 2.28e-12 ***
## AO1999.Mar         0.37531    0.09644   3.892 9.96e-05 ***
## AO2002.Oct         0.49217    0.08564   5.747 9.10e-09 ***
## AR-Nonseasonal-01  0.64477    0.07100   9.081  < 2e-16 ***
## MA-Seasonal-12     0.99965    0.08975  11.138  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 0 0)(0 1 1)  Obs.: 97  Transform: none
## AICc: -109.8, BIC: -88.36  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 20.02   Shapiro (normality): 0.9765 .

sliding <- series(m, "sfs")

## specs have been added to the model: slidingspans

sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100

## Warning in is.na(object): is.na() aplicado a um objeto diferente de lista
## ou vetor de tipo 'NULL'

## [1] NaN

#view(m)

capital_trabalho_2a <- final(m)

ggmonthplot(capital_trabalho_2a) + 
  theme_bw() +
  scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
  labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "De Jan/96 a Jan/04")

# Gráfico por ano e mês
ggseasonplot(capital_trabalho_2a, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
  labs(title = "Média mensal da razão capital-trabalho", x = "Mês", subtitle = "De Jan/96 a Jan/04")

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho_a), by = "1 month"), capital_trabalho_a, capital_trabalho_2a)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho_a, colour = "Original")) +
  geom_line(aes(y = capital_trabalho_2a, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Capital/Trabalho")+
  labs(color="Variável") +
  theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1

# Dividindo a série em pedaços distintos: 
# pega os pontos de quebras
bp_ts <- breakpoints(capital_trabalho_a ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = capital_trabalho_a ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1   24            
## m = 2   24 42         
## m = 3   24    49 65   
## m = 4   24    49 65 81
## m = 5   24 39 53 67 81
## 
## Corresponding to breakdates:
##                                                 
## m = 1   1997(12)                                
## m = 2   1997(12) 1999(6)                        
## m = 3   1997(12)         2000(1) 2001(5)        
## m = 4   1997(12)         2000(1) 2001(5) 2002(9)
## m = 5   1997(12) 1999(3) 2000(5) 2001(7) 2002(9)
## 
## Fit:
##                                              
## m   0      1      2      3      4      5     
## RSS  6.505  4.558  4.078  3.834  3.791  3.990
## BIC 22.318 -3.045 -4.687 -1.523  6.530 20.643

# Calcula o intervalo de confiança
ci_ts <- confint(bp_ts)


#############
m <- seas(x = capital_trabalho_b, transform.function = "none", regression.aictest = NULL)
qs(m)

##                     qs   p-val
## qsori        229.41027 0.00000
## qsorievadj   277.61642 0.00000
## qsrsd          0.00000 1.00000
## qssadj         3.02991 0.21982
## qssadjevadj    0.00000 1.00000
## qsirr          2.37126 0.30555
## qsirrevadj     0.00000 1.00000
## qssori       124.60633 0.00000
## qssorievadj  149.00319 0.00000
## qssrsd         0.00000 1.00000
## qsssadj        6.03848 0.04884
## qsssadjevadj   0.00000 1.00000
## qssirr         7.36388 0.02517
## qssirrevadj    0.00000 1.00000

plot(m)

summary(m)

## 
## Call:
## seas(x = capital_trabalho_b, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AO2004.Jun        -0.75270    0.08042  -9.359  < 2e-16 ***
## LS2004.Nov        -0.33098    0.04195  -7.889 3.03e-15 ***
## AO2004.Dec        -0.39821    0.05723  -6.959 3.44e-12 ***
## LS2005.Mar         0.15205    0.03782   4.020 5.82e-05 ***
## AO2005.Jun         0.37303    0.07240   5.152 2.57e-07 ***
## AO2006.Jun         0.51048    0.06305   8.097 5.63e-16 ***
## AO2006.Dec        -0.22518    0.04357  -5.168 2.36e-07 ***
## AO2007.Jun         0.35674    0.05244   6.803 1.02e-11 ***
## AO2008.Jan        -0.19683    0.04278  -4.601 4.21e-06 ***
## AO2010.Jun        -0.19763    0.04285  -4.612 3.98e-06 ***
## AO2011.Oct         0.23059    0.04291   5.374 7.69e-08 ***
## AO2011.Dec         0.22011    0.04291   5.130 2.90e-07 ***
## AO2013.Jun        -0.40288    0.04538  -8.878  < 2e-16 ***
## AO2014.Jun        -0.31472    0.04536  -6.939 3.96e-12 ***
## AO2016.Dec         0.20542    0.05233   3.925 8.67e-05 ***
## MA-Nonseasonal-01  0.73304    0.05551  13.204  < 2e-16 ***
## MA-Seasonal-12     0.33591    0.07718   4.352 1.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 1 1)(0 1 1)  Obs.: 166  Transform: none
## AICc: -403.3, BIC: -353.9  QS (no seasonality in final): 3.03  
## Box-Ljung (no autocorr.): 12.88   Shapiro (normality): 0.9618 ***

sliding <- series(m, "sfs")

## specs have been added to the model: slidingspans

sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100

## [1] 0

#view(m)
capital_trabalho_2b <- final(m)

ggmonthplot(capital_trabalho_2b) + 
  theme_bw() +
  scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
  labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "De Fev/04 a Nov/17")

# Gráfico por ano e mês
ggseasonplot(capital_trabalho_2b, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
  labs(title = "Média mensal da razão capital-trabalho", x = "Mês", subtitle = "De Fev/04 a Nov/17")

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("2004-02-01"), length = length(capital_trabalho_b), by = "1 month"), capital_trabalho_b, capital_trabalho_2b)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho_b, colour = "Original")) +
  geom_line(aes(y = capital_trabalho_2b, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Capital/Trabalho")+
  labs(color="Variável") +
  theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1

# Dividindo a série em pedaços distintos: 
# pega os pontos de quebras
bp_ts <- breakpoints(capital_trabalho_b ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = capital_trabalho_b ~ 1)
## 
## Breakpoints at observation number:
##                         
## m = 1                136
## m = 2   29           136
## m = 3   29 65        130
## m = 4   29 60 88     136
## m = 5   29 60 88 112 136
## 
## Corresponding to breakdates:
##                                                 
## m = 1                                   2015(5) 
## m = 2   2006(6)                         2015(5) 
## m = 3   2006(6) 2009(6)                 2014(11)
## m = 4   2006(6) 2009(1) 2011(5)         2015(5) 
## m = 5   2006(6) 2009(1) 2011(5) 2013(5) 2015(5) 
## 
## Fit:
##                                              
## m   0      1      2      3      4      5     
## RSS  13.57  13.33  13.12  13.09  13.02  13.02
## BIC  65.60  72.89  80.47  90.37  99.64 109.83

# Calcula o intervalo de confiança
#ci_ts <- confint(bp_ts)

###############
# Juntando novamente as duas séries
###############

comb <- ts.union(capital_trabalho_2a, capital_trabalho_2b)
capital_trabalho_final <- pmin(comb[,1], comb[,2], na.rm = TRUE)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho_final), by = "1 month"), capital_trabalho, capital_trabalho_final)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho, colour = "Original")) +
  geom_line(aes(y = capital_trabalho_final, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Razão Capital/Trabalho")+
  labs(color="Variável") +
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1

df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho2), by = "1 month"), capital_trabalho,  capital_trabalho2, capital_trabalho_final)
names(df) <- c("Data", "capital_trabalho", "capital_trabalho2", "capital_trabalho_final")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = capital_trabalho, colour = "Original")) +
  geom_line(aes(y = capital_trabalho2, colour = "Automático")) +
  geom_line(aes(y = capital_trabalho_final, colour = "Manual")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Razão Capital/Trabalho")+
  labs(color="Método") +
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1

#######################
# Selic               #
#######################
cores <- brewer.pal(6, "Dark2")
selic_4390 <- BETS.get("4390", from = "1996-01-01", to = "2017-11-30") 

## Transformando em série anual
selic <- ((1+selic_4390/100)^(12)-1)*100

ggmonthplot(selic) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a - média)") +
  labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")

# Gráfico por ano e mês
ggseasonplot(selic, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")

m <- seas(x = selic, transform.function = "none", regression.aictest = NULL) 
#plot(m)
summary(m)

## 
## Call:
## seas(x = selic, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## LS1996.Mar        -4.09295    0.97485  -4.199 2.69e-05 ***
## LS1997.Nov        21.72884    0.95722  22.700  < 2e-16 ***
## LS1998.Jan        -6.54038    0.99031  -6.604 3.99e-11 ***
## LS1998.Feb        -5.49776    0.98955  -5.556 2.76e-08 ***
## LS1998.Apr        -7.40308    0.95723  -7.734 1.04e-14 ***
## LS1998.Sep        15.32547    1.03446  14.815  < 2e-16 ***
## AO1998.Oct         7.50672    0.97560   7.694 1.42e-14 ***
## AO1998.Nov         3.92048    0.85822   4.568 4.92e-06 ***
## AO1999.Jan        -4.09824    0.92123  -4.449 8.64e-06 ***
## AO1999.Mar        14.93677    0.85625  17.444  < 2e-16 ***
## LS1999.May        -5.36681    1.05880  -5.069 4.00e-07 ***
## LS1999.Jun        -4.48989    0.92662  -4.845 1.26e-06 ***
## AR-Nonseasonal-01 -0.25889    0.10018  -2.584  0.00976 ** 
## AR-Nonseasonal-02  0.13521    0.06290   2.150  0.03159 *  
## AR-Nonseasonal-03  0.49771    0.05109   9.742  < 2e-16 ***
## MA-Nonseasonal-01  0.14349    0.11460   1.252  0.21053    
## MA-Seasonal-12     1.00000    0.04180  23.923  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 263  Transform: none
## AICc: 799.2, BIC: 859.7  QS (no seasonality in final):2.135  
## Box-Ljung (no autocorr.): 59.74 *** Shapiro (normality): 0.9856 **
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.
## 
## Notes:
## - Unable to test LS1998.Jan due to regression matrix
##   singularity.
## - Unable to test LS1999.Jun due to regression matrix
##   singularity.
## - Unable to test LS1998.Feb due to regression matrix
##   singularity.
## - Unable to test LS1998.Feb due to regression matrix
##   singularity.
## - Unable to test AO1998.Jan due to regression matrix
##   singularity.
## - Unable to test LS1999.May due to regression matrix
##   singularity.

qs(m)

##                    qs   p-val
## qsori        23.76401 0.00001
## qsorievadj   96.16089 0.00000
## qsrsd         0.00000 1.00000
## qssadj        2.13450 0.34395
## qssadjevadj  18.30242 0.00011
## qsirr         0.00000 1.00000
## qsirrevadj    0.00000 1.00000
## qssori       33.30455 0.00000
## qssorievadj  33.30455 0.00000
## qssrsd        0.26400 0.87634
## qsssadj       6.61682 0.03657
## qsssadjevadj  6.61682 0.03657
## qssirr        0.00000 1.00000
## qssirrevadj   0.00000 1.00000

#summary(m) diz Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks se for a série a partir de 99

# Usa a interface gráfica
#view(m)
sliding <- series(m, "sfs")

## specs have been added to the model: slidingspans

sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100

## [1] 0

# Salva a série final em uma nova variável
selic2 <- final(m)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date(inicio), length = length(selic2), by = "1 month"), selic, selic2)
names(df) <- c("Data", "Original", "Série Filtrada")
p2 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = selic, colour = "Original")) +
  geom_line(aes(y = selic2, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Selic (%a.a.)")+
  labs(color="Variável") +
  theme_bw()+theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p2

# Dividindo a série em pedaços distintos: 
# pega os pontos de quebras
bp_ts <- breakpoints(selic ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = selic ~ 1)
## 
## Breakpoints at observation number:
##                          
## m = 1         121        
## m = 2   41    128        
## m = 3   41    128     222
## m = 4   41    123 162 222
## m = 5   41 81 121 160 222
## 
## Corresponding to breakdates:
##                                                
## m = 1                   2006(1)                
## m = 2   1999(5)         2006(8)                
## m = 3   1999(5)         2006(8)         2014(6)
## m = 4   1999(5)         2006(3) 2009(6) 2014(6)
## m = 5   1999(5) 2002(9) 2006(1) 2009(4) 2014(6)
## 
## Fit:
##                                        
## m   0     1     2     3     4     5    
## RSS 13285  5841  3529  3417  3262  3211
## BIC  1789  1584  1463  1465  1464  1471

# Calcula o intervalo de confiança
ci_ts <- confint(bp_ts)

Data <- seq(as.Date("1996-01-01"), length = length(selic), by = "1 month")
df1 <- melt(data.frame(Data,selic), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")

p1 <- ggplot(df1[which(df1$Variavel == "selic"),], aes(Data, Valor, colour = Variavel)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Selic (%a.a.)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+ 
        geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
        geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(selic), xend = Data[ci_ts$confint[5]], yend = min(selic)))+
        geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(selic), xend = Data[ci_ts$confint[6]], yend = min(selic)))+
        theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1

###########
# Dividindo a selic em três pedaços
##########

# Salvando três séries para passar o seasonal
selic_a <- selic[1:(bp_ts$breakpoints[1]-1)]
selic_b <- selic[(bp_ts$breakpoints[1]):(bp_ts$breakpoints[2]-1)]
selic_c <- selic[(bp_ts$breakpoints[2]):(length(selic))]
selic_a <- ts(selic_a,  start = c(1996, 01, 01), frequency = 12)
selic_b <- ts(selic_b,  start = c(1999, 05, 01), frequency = 12)
selic_c <- ts(selic_c,  start = c(2006, 08, 01), frequency = 12)

#######
# Primeiro pedaço
#########

ggmonthplot(selic_a) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a - média)") +
  labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Jan/96 a Abr/99")

# Gráfico por ano e mês
ggseasonplot(selic_a, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Jan/96 a Abr/99")

m <- seas(x = selic_a, transform.function = "none", regression.aictest = NULL) 
plot(m)

summary(m)

## 
## Call:
## seas(x = selic_a, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## Constant        -0.6762     0.3314  -2.040   0.0413 *  
## LS1997.Nov      21.4637     1.3511  15.886  < 2e-16 ***
## LS1998.Feb      -9.0236     1.2976  -6.954 3.55e-12 ***
## LS1998.Apr      -7.0057     1.3511  -5.185 2.16e-07 ***
## LS1998.Sep      15.4392     1.4204  10.870  < 2e-16 ***
## AO1998.Oct       5.8520     1.0364   5.646 1.64e-08 ***
## AO1999.Mar      15.4913     1.0364  14.947  < 2e-16 ***
## AR-Seasonal-12   0.4275     0.1998   2.140   0.0323 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 1 0)(1 0 0)  Obs.: 40  Transform: none
## AICc:   164, BIC: 172.7  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 31.97   Shapiro (normality): 0.9405 *
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the prior adjusted series (Table
##   B1) has no visually significant seasonal peaks.
## 
## Notes:
## - Insufficient data to compute average forecast error
##   diagnostic.
## - The program cannot perform hypothesis tests for kurtosis
##   on less than 50 observations.

qs(m)

##                  qs   p-val
## qsori       0.02853 0.98584
## qsorievadj  4.23955 0.12006
## qsrsd       1.32801 0.51479
## qssadj      0.00000 1.00000
## qssadjevadj 0.00000 1.00000
## qsirr       0.01122 0.99441
## qsirrevadj  0.00000 1.00000

#summary(m) #diz Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks

# Deixei ela como está

selic_2a <- selic_a

#######
# Segundo pedaço
#########

ggmonthplot(selic_b) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a - média)") +
  labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")

# Gráfico por ano e mês
ggseasonplot(selic_b, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")

m <- seas(x = selic_b, transform.function = "none", regression.aictest = NULL) 
plot(m)

summary(m)

## 
## Call:
## seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -0.19030    0.14589  -1.304    0.192    
## AR-Nonseasonal-02  0.09671    0.09194   1.052    0.293    
## AR-Nonseasonal-03  0.51926    0.08351   6.218 5.03e-10 ***
## MA-Nonseasonal-01  0.12082    0.17475   0.691    0.489    
## MA-Seasonal-12     0.99913    0.12930   7.727 1.10e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 87  Transform: none
## AICc: 306.6, BIC: 319.2  QS (no seasonality in final):3.007  
## Box-Ljung (no autocorr.): 44.68 ** Shapiro (normality): 0.9865  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.

qs(m)

##                   qs   p-val
## qsori       30.28836 0.00000
## qsorievadj  30.28836 0.00000
## qsrsd        0.00000 1.00000
## qssadj       3.00726 0.22232
## qssadjevadj  3.00726 0.22232
## qsirr        0.00000 1.00000
## qsirrevadj   0.00000 1.00000

summary(m)

## 
## Call:
## seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -0.19030    0.14589  -1.304    0.192    
## AR-Nonseasonal-02  0.09671    0.09194   1.052    0.293    
## AR-Nonseasonal-03  0.51926    0.08351   6.218 5.03e-10 ***
## MA-Nonseasonal-01  0.12082    0.17475   0.691    0.489    
## MA-Seasonal-12     0.99913    0.12930   7.727 1.10e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 87  Transform: none
## AICc: 306.6, BIC: 319.2  QS (no seasonality in final):3.007  
## Box-Ljung (no autocorr.): 44.68 ** Shapiro (normality): 0.9865  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.

selic_2b <- final(m)

ggmonthplot(selic_2b) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média diária e mensal da Taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")

# Gráfico por ano e mês
ggseasonplot(selic_2b, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média mensal da Taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")

selic_2b <- final(m)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1999-05-01"), length = length(selic_b), by = "1 month"), selic_b, selic_2b)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = selic_b, colour = "Original")) +
  geom_line(aes(y = selic_2b, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Selic (%a.a.)")+
  labs(color="Série") +
  theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1

#############
# Terceiro pedaço
#############

ggmonthplot(selic_c) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a - média)") +
  labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")

# Gráfico por ano e mês
ggseasonplot(selic_c, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")

m <- seas(x = selic_c, transform.function = "none", regression.aictest = NULL) 
plot(m)

summary(m)

## 
## Call:
## seas(x = selic_c, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -0.15177    0.10739  -1.413 0.157587    
## AR-Nonseasonal-02  0.23574    0.07095   3.323 0.000891 ***
## AR-Nonseasonal-03  0.58420    0.06552   8.916  < 2e-16 ***
## MA-Nonseasonal-01  0.31636    0.12889   2.455 0.014107 *  
## MA-Seasonal-12     0.99998    0.08835  11.319  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 136  Transform: none
## AICc: 305.7, BIC: 321.9  QS (no seasonality in final): 2.96  
## Box-Ljung (no autocorr.): 57.69 *** Shapiro (normality): 0.985  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.

qs(m)

##                    qs   p-val
## qsori        44.66623 0.00000
## qsorievadj   44.66623 0.00000
## qsrsd         0.00000 1.00000
## qssadj        2.96000 0.22764
## qssadjevadj   2.96000 0.22764
## qsirr         0.00000 1.00000
## qsirrevadj    0.00000 1.00000
## qssori       33.30455 0.00000
## qssorievadj  33.30455 0.00000
## qssrsd        0.00000 1.00000
## qsssadj       3.39334 0.18329
## qsssadjevadj  3.39334 0.18329
## qssirr        0.00000 1.00000
## qssirrevadj   0.00000 1.00000

selic_2c <- final(m)

ggmonthplot(selic_2c) + 
  theme_bw() +
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média diária e mensal da Taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")

# Gráfico por ano e mês
ggseasonplot(selic_2c, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Selic (%a.a. - média)") +
  labs(title = "Média mensal da Taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("2006-08-01"), length = length(selic_c), by = "1 month"), selic_c, selic_2c)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = selic_c, colour = "Original")) +
  geom_line(aes(y = selic_2c, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Selic (%a.a.)")+
  labs(color="Série") +
  theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1

###############
# Juntando novamente as três séries
###############

comb <- ts.union(selic_2a, selic_2b, selic_2c)
selic_final <- pmin(comb[,1], comb[,2], comb[,3], na.rm = TRUE)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(selic_final), by = "1 month"), selic, selic_final)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = selic, colour = "Original")) +
  geom_line(aes(y = selic_final, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Selic (%a.a.)")+
  labs(color="Série") +
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1

m <- seas(x = selic, transform.function = "none", regression.aictest = NULL)
summary(m)

## 
## Call:
## seas(x = selic, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## LS1996.Mar        -4.09295    0.97485  -4.199 2.69e-05 ***
## LS1997.Nov        21.72884    0.95722  22.700  < 2e-16 ***
## LS1998.Jan        -6.54038    0.99031  -6.604 3.99e-11 ***
## LS1998.Feb        -5.49776    0.98955  -5.556 2.76e-08 ***
## LS1998.Apr        -7.40308    0.95723  -7.734 1.04e-14 ***
## LS1998.Sep        15.32547    1.03446  14.815  < 2e-16 ***
## AO1998.Oct         7.50672    0.97560   7.694 1.42e-14 ***
## AO1998.Nov         3.92048    0.85822   4.568 4.92e-06 ***
## AO1999.Jan        -4.09824    0.92123  -4.449 8.64e-06 ***
## AO1999.Mar        14.93677    0.85625  17.444  < 2e-16 ***
## LS1999.May        -5.36681    1.05880  -5.069 4.00e-07 ***
## LS1999.Jun        -4.48989    0.92662  -4.845 1.26e-06 ***
## AR-Nonseasonal-01 -0.25889    0.10018  -2.584  0.00976 ** 
## AR-Nonseasonal-02  0.13521    0.06290   2.150  0.03159 *  
## AR-Nonseasonal-03  0.49771    0.05109   9.742  < 2e-16 ***
## MA-Nonseasonal-01  0.14349    0.11460   1.252  0.21053    
## MA-Seasonal-12     1.00000    0.04180  23.923  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 1 1)(0 1 1)  Obs.: 263  Transform: none
## AICc: 799.2, BIC: 859.7  QS (no seasonality in final):2.135  
## Box-Ljung (no autocorr.): 59.74 *** Shapiro (normality): 0.9856 **
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.
## 
## Notes:
## - Unable to test LS1998.Jan due to regression matrix
##   singularity.
## - Unable to test LS1999.Jun due to regression matrix
##   singularity.
## - Unable to test LS1998.Feb due to regression matrix
##   singularity.
## - Unable to test LS1998.Feb due to regression matrix
##   singularity.
## - Unable to test AO1998.Jan due to regression matrix
##   singularity.
## - Unable to test LS1999.May due to regression matrix
##   singularity.

selic2 <- final(m)

df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho2), by = "1 month"), selic,  selic2, selic_final)
names(df) <- c("Data", "selic", "selic2", "selic_final")
p1 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = selic, colour = "Original")) +
  geom_line(aes(y = selic2, colour = "Automático")) +
  geom_line(aes(y = selic_final, colour = "Manual")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "Selic (%a.a.)")+
  labs(color="Método") +
  theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1

#######################
# IBC-Br              #
#######################

m <- seas(x = ibcbr)
#plot(m)
#summary(m)

# Usa a interface gráfica
#view(m)

# Como ficou igual, também não mexi nele.

#######################
# PIB Mensal          #
#######################
cores <- brewer.pal(6, "Dark2")
ggmonthplot(PIB) + 
  theme_bw() +
  scale_y_continuous(name = "PIB (log-dif - média)") +
  labs(title = "Média diária e mensal da variação do PIB", x = "Mês", subtitle = "Nenhuma transformação")

# Gráfico por ano e mês
ggseasonplot(PIB, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "PIB (log-dif - média)") +
  labs(title = "Média mensal da da variação do PIB", x = "Mês", subtitle = "Nenhuma transformação")

m <- seas(x = PIB, transform.function = "none", regression.aictest = NULL)

## Model used in SEATS is different: (2 0 2)(0 1 1)

plot(m)

summary(m)

## 
## Call:
## seas(x = PIB, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AR-Nonseasonal-01 -1.10213    0.06988 -15.771  < 2e-16 ***
## AR-Nonseasonal-02 -0.38179    0.08846  -4.316 1.59e-05 ***
## AR-Nonseasonal-03 -0.13937    0.06474  -2.153   0.0313 *  
## MA-Nonseasonal-01 -0.96137    0.03241 -29.666  < 2e-16 ***
## MA-Seasonal-12     0.79111    0.04043  19.568  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (3 0 1)(0 1 1)  Obs.: 263  Transform: none
## AICc: -1192, BIC: -1171  QS (no seasonality in final):0.8726  
## Box-Ljung (no autocorr.): 86.12 *** Shapiro (normality): 0.9961  
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
##   been found in one or more of the estimated spectra.
## 
## Notes:
## - Model used for SEATS decomposition is different from the
##   model estimated in the regARIMA modeling module of
##   X-13ARIMA-SEATS.

# Usa a interface gráfica
#view(m)

# Salva a série final em uma nova variável
PIB2 <- final(m)

# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(PIB), by = "1 month"), PIB, PIB2)
names(df) <- c("Data", "Original", "Série Filtrada")
p2 <- ggplot(df, aes(Data)) +
  geom_line(aes(y = PIB, colour = "Original")) +
  geom_line(aes(y = PIB2, colour = "Filtrada")) + 
  scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
  scale_y_continuous(name = "PIB Mensal (variação mensal do índice, base = 11/17)")+
  labs(color="Variável") +
  theme_bw()
p2 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 8))

#######################
# Câmbio              #
#######################
cores <- brewer.pal(6, "Dark2")
ggmonthplot(cambio) + 
  theme_bw() +
  scale_y_continuous(name = "Tx. Cambio (log-dif - média)") +
  labs(title = "Média diária e mensal da variação do Câmbio", x = "Mês", subtitle = "Nenhuma transformação")

# Gráfico por ano e mês
ggseasonplot(PIB, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "Tx. Cambio (log-dif - média)") +
  labs(title = "Média mensal da variação do Câmbio", x = "Mês", subtitle = "Nenhuma transformação")

m <- seas(x = cambio, transform.function = "none", regression.aictest = NULL)
plot(m)

summary(m) # Summary diz que não deve fazer ajuste sazonal

## 
## Call:
## seas(x = cambio, transform.function = "none", regression.aictest = NULL)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## AO1999.Jan         0.48976    0.04001  12.241  < 2e-16 ***
## AO1999.Mar        -0.18202    0.04001  -4.549 5.38e-06 ***
## AO2002.Jul         0.18745    0.04001   4.685 2.80e-06 ***
## AO2002.Sep         0.27601    0.04001   6.899 5.25e-12 ***
## AO2011.Sep         0.16420    0.04001   4.104 4.06e-05 ***
## AR-Nonseasonal-01  0.11781    0.06123   1.924   0.0544 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (1 0 0)  Obs.: 263  Transform: none
## AICc: -928.6, BIC:  -904  QS (no seasonality in final):0.4173  
## Box-Ljung (no autocorr.): 31.84   Shapiro (normality): 0.9713 ***
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the prior adjusted series (Table
##   B1) has no visually significant seasonal peaks.

# Usa a interface gráfica
#view(m)

#######################
# IPCA                #
#######################
cores <- brewer.pal(6, "Dark2")
ggmonthplot(ipca) + 
  theme_bw() +
  scale_y_continuous(name = "IPCA (% Acum. 12 meses - média)") +
  labs(title = "Média diária e mensal do IPCA acum. 12 meses", x = "Mês", subtitle = "Nenhuma transformação")

# Gráfico por ano e mês
ggseasonplot(ipca, year.labels = TRUE) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(name = "IPCA (% Acum. 12 meses - média)") +
  labs(title = "Média mensal do IPCA acum. 12 meses", x = "Mês", subtitle = "Nenhuma transformação")

m <- seas(x = ipca, transform.function = "none", regression.aictest = NULL)

## Model used in SEATS is different: (0 2 1)

#plot(m)
#summary(m) # Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks.

# Usa a interface gráfica
#view(m)

# Repete o mesmo gráfico de antes, mas com as variáveis dessazonalizadas

df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho2, selic2, ibcbr, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "IBCBr", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")

# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Capital trabalho") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p2 <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[2])+
        scale_y_continuous(name="Selic (%a.a.)") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p3 <- ggplot(df2[which(df2$variable == "IBCBr"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[3])+
        scale_y_continuous(name="IBC-Br") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[4])+
        scale_y_continuous(name="Tx. Câmbio (var)") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p5 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[5])+
        scale_y_continuous(name="IPCA (acum. 12m.)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
        theme_bw()
p5 <- p5 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5)

# Repete o mesmo gráfico de antes, mas com as variáveis dessazonalizadas
cores <- brewer.pal(6, "Dark2")
df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho_final, selic_final, PIB2, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "PIB", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")

## Warning: attributes are not identical across measure variables; they will
## be dropped

# Gráficos individuais
p1a <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Capital trabalho") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p1a <- p1a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p2a <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[2])+
        scale_y_continuous(name="Selic (%a.a.)") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p2a <- p2a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p3a <- ggplot(df2[which(df2$variable == "PIB"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[3])+
        scale_y_continuous(name="PIB") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p3a <- p3a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p4a <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[4])+
        scale_y_continuous(name="Tx. Câmbio (var)") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p4a <- p4a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p5a <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[5])+
        scale_y_continuous(name="IPCA (acum. 12m.)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
        theme_bw()
p5a <- p5a + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

grid.arrange(p1a, p2a, p3a, p4a, p5a, ncol=1, nrow = 5)

Parece que o X-13ARIMA fez um trabalho melhorzinho, então vou continuar com a série capital_trabalho_adj.

Descritivas

Fiz as descritivas para a série original e depois para a série trabalhada manualmente e pelo X-13ARIMA. Fiz as descritivas das variáveis que vão entrar no VAR. Depois tem que lembrar de fazer das que efetivamente ficaram Uma ideia é colocar naqueles arquivos dos VAR individuais. Acho que vou fazer isso depois…

### Descriptives

descriptives     <- matrix(NA, nrow = 8, ncol = (ncol(df1)-1))
rownames(descriptives) <- c("Observações", "Mínimo", "1o quartil",
                      "Média", "Mediana",  "3o quartil", "Máximo",
                      "Desv. Pad.")

colnames(descriptives) <- names(df1)[-1]

desc <- function(x) {
  n       <- length(x)
  minimum <- min(x, na.rm = TRUE)
  first_q <- quantile(x, 0.25, na.rm = TRUE)
  media   <- mean(x, na.rm = TRUE)
  mediana <- median(x, na.rm = TRUE)
  third_q <- quantile(x, 0.75, na.rm = TRUE)
  maximum <- max(x, na.rm = TRUE)
  std     <- sd(x, na.rm = TRUE)

    return(list(n = n, minimum = minimum, first_quar = first_q, media = media, mediana = mediana, third_quar = third_q, maximum = maximum, std = std))
}

for (i in 1:8){
  descriptives[i, 1] <- round(as.numeric(desc(df1[,2])[i]),4)
  descriptives[i, 2] <- round(as.numeric(desc(df1[,3])[i]),4)
  descriptives[i, 3] <- round(as.numeric(desc(df1[,4])[i]),4)
  descriptives[i, 4] <- round(as.numeric(desc(df1[,5])[i]),4)
  descriptives[i, 5] <- round(as.numeric(desc(df1[,6])[i]),4)
}

descriptives[1,] <- as.integer(descriptives[1,])

descriptives <- data.frame(descriptives)

#stargazer(descriptives, summary=FALSE, header = TRUE, type = 'html')
stargazer(descriptives, summary=FALSE, header = TRUE, type = 'latex')

% Table created by stargazer v.5.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu % Date and time: qua, jan 24, 2018 - 03:16:31

Baixando dados de desigualdade

Lista de dados existentes

Jogando no liquidificador

O pacote que tem o modelo implementado é o bvarsv. Aqui (https://github.com/FK83/bvarsv/blob/master/bvarsv_Nov2015_website.pdf) e aqui (https://github.com/FK83/bvarsv/blob/master/bvarsv_replication.pdf) tem exemplos de uso.

Especificações da função bvar.sv.tvp:

• \(p\) é o número de defasagens;
• \(\tau\) é o número de observações que serão utilizadas para conseguir os hiperparâmetros. O default é 40, usei 24 (2 anos);
• \(nf\) é o número de períodos para fazer previsão, o default é 10;
• $pdrift* parâmetro de tendência, default é TRUE;
• As distribuições a priori são as mesmas do Primiceri:
  • \(B_0\) são os betas iniciais e seguem uma distribuição normal com parâmetros obtidos usando MQO
  • \(A_0\) é a covariância inicial, também tem distribuição normal. As variâncias de \(A_0\) e \(B_0\) são multiplicadas por 4;
  • \(log \ \sigma_0\) é a log volatilidade inicial, também segue uma densidade normal, com média estimada por MQO e variância unitária;
  • \(Q\) matriz de variância covariância de \(B_t\), tem distribuição inversa-wishart e uns parâmetros meio esquisitos (https://github.com/FK83/bvarsv/blob/master/bvarsv_Nov2015_website.pdf - ver página 2)
  • \(W\) é a matriz decovariância dos choques em \(log\ \sigma_0\)
  • \(S_j\) é a matriz de covariância de \(A_t\)

Coisas a serem pensadas

• Uma coisa a ser feita é jogar esse IBC-Br fora e aumentar um pouco a amostra para ver o que acontece.
• Tem a questão do câmbio também pra ser vista.
• E a questão da priori, talvez mude os resultados se eu mudar ela.

Referências